Debye-Hiickel-Bjerrum theory for charged colloids 
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We formulate an extension of the Debye-Hiickel-Bjerrum theory [M. E. Fisher and Y. Levin, Phys. 
Rev. Lett. 71, 3826 (1993)] to the fluid state of a highly asymmetric charged colloid. Allowing for 
the formation of clusters consisting of one polyion and n condensed counterions, the total Helmholtz 
free energy of the colloidal suspension is constructed. The thermodynamic properties, such as the 
cluster-density distribution and the pressure, are obtained by the minimization of the free energy 
under the constraints of fixed number of polyions and counterions. In agreement with the current 
experimental and Monte Carlo results, no evidence of any phase transition is encountered. 

PACS numbers: 82.70.Dd; 36.20.-r; 64.60.Cn 



I. INTRODUCTION 

The technological importance of charged colloidal sus- 
pensions can not be overemphasized. One comes across 
these important systems in fields as diverse as the chem- 
ical engineering and the environmental science. Many 
water-soluble paints contain charged colloidal suspen- 
sions as a main ingredient. In this case the problem 
of great industrial importance is to stabilize the suspen- 
sion against the flocculation and precipitation. On the 
other extreme is the constantly growing environmental 
necessity of cleaning contaminated water. For this it is 
essential to find the most effective way of precipitating 
the (usually) charged organic particles dissolved in the 
water. 

From the theoretical perspective the problem of 
strongly asymmetric electrolyte solutions is extremely 
difficult to study. The long-ranged nature of the Coulomb 
force combined with the large charge and size asymme- 
try between the polyions and the counterions and coions 
makes it impossible to use the traditional methods of 
liquid-state theory. At high volume fractions the suspen- 
sion will crystallize, that is, the polyions will become ar- 
ranged in the form of a lattice. The solid state provides a 
major simplification of reducing the many-polyion prob- 
lem to that of one polyion inside a Wigner-Seitz (WS) cell 
M. Unfortunately, at low densities, or in the presence of 
a simple electrolyte, the suspension becomes disordered 
and the WS picture is no longer valid . A new strategy 
must be tried. 

The Debye-Hiickel-Bjerrum theory (DHBj) || was 
quite successful in explaining the behavior of symmet- 
ric electrolytes. The Bjerrum's concept of association of 
oppositely charged ions into dipolar pairs || served to 
correct the Debye-Hiickel (DH) linearization of the non- 
linear Poisson-Boltzmann equation pi. Taking into ac- 



count the dipolar solvation energy made the coexistence 
curve produced by the DHBj theory ||] become in excel- 
lent agreement with the Monte Carlo simulations || . The 
large surface charge of a polyion suggests that the cluster 
formation should be even more important in the case of 
polyelectrolytes. In the present work we shall explore to 
what extent the counterion condensation influences the 
thermodynamic properties of a polyelectrolyte solution. 
The DHBj theory will be extended to treat the fluid state 
of a charged colloidal suspension. 



II. DEFINITION OF THE MODEL 



We shall work with the primitive model of polyelec- 
trolyte (PMP) jjj. The system will consist of N p polyions 
inside a volume V. The polyions are modeled as hard 
spheres of radius a with a uniform surface charge density 
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where Z is the polyion valence (number of ionized sites) 
and q > is the charge of a proton. To preserve the 
overall charge neutrality of the system, ZN P point-like 
counterions of charge +q are present. The solvent is mod- 
eled as a homogeneous medium of dielectric constant D. 
Due to the strong electrostatic interaction between the 
polyions and the counterions, we expect that the asym- 
metric polyelectrolyte will be composed of clusters, with 
density p n , consisting of one polyion and 1 < n < Z asso- 
ciated counterions, as well as bare polyions of density po 
and free (unassociated) counterions of density pf. The 
conservation of the total number of polyions and counte- 
rions leads to two conservation equations, 
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n=0 



and 



(3) 



where p p — N p /V is the total density of polyions (asso- 
ciated or not). 

All the thermodynamic properties of the system can 
be determined once the free energy is known. The 
Helmholtz free energy can be split into two parts, the 
electrostatic and entropic. The electrostatic terms are 
due to the inter-particles interactions and can be at- 
tributed to the polyion-counterion, the polyion-polyion, 
and the counterion-counterion interactions. All the elec- 
trostatic interactions will be evaluated using the DH the- 
ory H . This is motivated by the former success of the 
theory when it was applied to symmetric electrolytes. In 
principle, the linear DH theory is satisfactory only for 
low densities and high temperatures. However, once the 
concept of clusters is introduced, the validity of the DHBj 
theory is extended into the nonlinear regime |]. The lin- 
ear structure of the DHBj theory insures its internal self 
consistency, a problem which is intrinsic to many of the 
nonlinear theories of electrolyte solutions ||. 

We shall assume that the effect of the counterion con- 
densation is to renormalize the polyion charge. Thus, the 
effective surface charge of a n-cluster is 
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All the nonlincarities related to the internal degrees of 
freedom of the clusters will be included in the entropic 
terms. In a previous work [7a] we have considered that 
the bounded counterions condense onto the surface of the 
spherical polyion. Although less realistic, this assump- 
tion has allowed us to obtain closed analytical expressions 
for the entropic contribution. In the present work, how- 
ever, the intra-cluster interactions will be treated using 
a local-density functional theory, so that the correlations 
between the bound counterions are explicitly taken into 
account. These correlations effects can be disregarded 
only when the concentration of counterions is not too 
high, a condition which may not be fulfilled in the close 
vicinity of a highly charged polyion. We now proceed to 
describe each one of the contributions to the Helmholtz 
free-energy density / = —F/V. 



III. THE POLYION-COUNTERION 
INTERACTION 

The polyion-counterion contribution is obtained using 
the usual DH theory applied to a n-cluster of effective 



surface charge a n inside the ionic atmosphere [7a]. Con- 
sider a n-cluster fixed at the origin, r — 0. Due to the 
hard-sphere exclusion, no counterions will be found in- 
side the region r < a, that is, 



p q (r < a) = . 



(5) 



Outside the spherical polyion, r > a, the cluster- 
counterion correlation function is approximated by a 
Boltzmann factor, leading to the charge density 



p g (r >a) = qp f exp [-0qip n (r)] -^(Z - n)qp n 

n=0 

+a n S(r-a), (6) 



where /3 _1 = /cbT and ip n (i") is the electrostatic poten- 
tial at a distance r from the center of the polyion. Notice 
that only the free counterions are assumed to get polar- 
ized; the bare polyions and clusters are too massive to be 
affected by the electrostatic fluctuations and contribute 
only to the neutralizing background. Substituting the 
charge density into the Poisson equation, 
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one obtains the nonlinear Poisson-Boltzmann equation. 
After the linearization of the exponential factor in (Q), 
the electrostatic potential VVi( r ) satisfies the Laplace (for 
r < a) and the Helmholtz (for r > a) equations, 



V 2 V n (r) 



, for r < a , 

K 2 ipn(r) - TJ^n <5( r - a) , for r > a , 



(8) 



where k = ^/47tAb/0/ is the inverse Debye screening 
length, and Ab = (3q 2 /D is the Bjerrum length. In prin- 
ciple the linearization is valid only in the limit pqip n <C 1, 
however, since the formation of clusters is taken into ac- 
count, the validity of the theory is extended into the non- 
linear regime ||. 

The second-order differential equation for ?Aii( r ) can 
be integrated, supplemented by the boundary conditions 
of vanishing of the electrostatic potential at infinity, the 
continuity of ipn{r) at r = a, and the discontinuity in 
the radial component of the electric field related to the 
presence of the surface charge a n at r — a. Under these 
conditions, we obtain 
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Using the charge density in the linearized form, 



Pq(r) = 



, for r < a , 

-^-ipnir) + o n S(r - a) , for r > a , 
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the electrostatic energy of a ro-cluster is calculated to be 



U n (K,q) = - I »AV/»,,(r)r„(rl 
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The electrostatic free-energy density for the polyion- 
counterion interaction is obtained through the Debye 
charging process, in which all the particles are simulta- 
neously charged from to their final charge p], 
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IV. THE POLYION-POLYION INTERACTION 



Due to the large asymmetry between the polyions 
and the counterions, the degrees of freedom associated 
with the counterions can be effectively integrated out. 
The long-ranged interaction between two clusters will be 
screened by the cloud of free counterions, producing an 
effective short-rang ed potential of a DLVO form g,[Kj], 



VZ(r)=q*(Z-n)(Z-m) 
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The polyion-polyion contribution to the Helmholtz free 
energy can then be calculated in the spirit of the usual 
van-der-Waals theory [3,7a], 
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V. THE COUNTERION-COUNTERION 
INTERACTION 



The counterion-counterion contribution, originating 
from the interactions between the free counterions, is 
calculated using the One Component Plasma (OCP) the- 
ory |Tl|j . The electrostatic free energy is found through 
a Debye charging process and a closed analytic form for 
f (pf), valid over a wide range of densities, is presented 
in [lib], 
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In the bulk this contribution is very small, and is included 
only for completeness. 



VI. THE MIXING FREE ENERGY 

The mixing free energy reduces to a sum of ideal-gas 
terms, 



Pf mbC (pf,iPn}) = E P> t 1 - ln (M'/C-)] , (18) 

where s 6 {f',n = Q,...,Z}, A s are the thermal de 
Broglie wavelengths associated with free counterions, 
bare polyions, and clusters; Cs are the internal parti- 
tion functions for an isolated specie s. Since the bare 
polyions and the free (unassociated) counterions do not 
have internal structure, their internal partition functions 
are simply given by £o = C/ = 1- F° r a n-cluster the 
internal partition function is 
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where the integration hypervolume fl n = {a < \n\ < 
R ni Vi = 1, • • • , n} depends on the cutoff R n . To fix the 
value of the cutoff, we follow an argument similar to the 
one used by Bjerrum in his study of dipolar formation in 
simple electrolytes 0. Suppose that we have a (n — 1)- 
cluster and we want to condense one more counterion 
to form a n-cluster. Because of the spherical symme- 
try, the (n—1) bound counterions contribute only to the 
renormalization of the polyion charge. The probability of 
finding the n th counterion at a distance r in the interval 
dr is 
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P{r) dr oc drr 2 exp [—0q<fi(r)] 

= dr r 2 exp \{Z — n + 1) Xb/t] , 



(21) 



where 4>{r) is the electrostatic potential generated by the 
(n — l)-cluster. The probability distribution P(r) has a 
minimum at 

A B 



r = R n = (Z-n+1)- 



(22) 



which, following Bjerrum M, we shall interpret as the 
distance of closest approach at which the n th counterion 
will become bound to the (n— l)-cluster. Since R n /a > 1, 
for a given reduced temperature, T* = a/Ae, there is a 
minimum value of the valence, Z min = IT* — 1, below 
which no counterions can condense onto a polyion, that 
is, the thermal entropic energy, 2k^T, will overcome the 
gain in electrostatic potential energy, (Z — n + l)q 2 /(Da), 
preventing the confinement from taking effect. 



With the cutoff defined, we shall now attempt to 
calculate the internal partition function of a n-cluster, 
Eq. (19). That, in itself, is a formidable task, since it 
requires evaluation of the many-body integrals (0). In- 
stead of performing the integrations explicitly, we shall 
use the local-density-functional theory to find the free 
energy of the condensed layer of n counterions, (3F£° n = 
— In C„ . Let us define the local density of counterions in 
the condensed layer of a n-cluster as 
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Within the local-density approximation (LDA), the 
Helmholtz free-energy functional f3J 7 ^ l ° n [g c (r)] corre- 
sponding to density g c (r) is 



(3f°°»[Q c (r)}= f d 3 rg c (r){ln[g c (r)Al] - 1} - ZA B f d 3 r 
Jv n Jv„ 

+ d 3 r g c {r)T co „[g c {r)] , 
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where the integrations are over the annulus V n = {a < 
\f \ < R n } and J- CO n is given by Eq. ([l6|). The first term 
in ( p4[ ) corresponds to the usual ideal-gas contribution, 
the second and the third terms are due to the electro- 
static interactions and the last term is the result of the 
correlations between the bounded counterions, for which 
we use the expression of the OCP theory |ll| . The equi- 
librium configuration, p c (r) = (g c {r)), is the one that 
minimizes the free-energy functional /3jF™ n [g c (r)] under 
the constraint 



d 3 r p c (r) = n . 



(25) 



This minimization procedure leads to the Boltzmann dis 
tribution for the density profile, 

nexp [-jicorrM - M{r)\ 



Pc{r) = 



d 3 r' exp [-Mcorr(r') - M{r')} 



(26) 



where the electrostatic and the excess chemical potentials 
are given, respectively, by 



orr (r) = T COTr [p c (r)] + p c (r) 
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with to given by Eq. ([IT]). On the other hand, the electro- 
static potential and the total charge density satisfy the 
Poisson equation, 



VVO) = V ■ E(r) 
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where E(r) — —Vil>(r) is the electric field at the point 
r. Inserting ( p6[ ) into (g9|) we find a Poisson-Boltzmann- 
like equation, which is a second-order nonlinear dif- 
ferential equation for the electrostatic potential ip(r). 
It should be remarked that, neglecting the correlation 
term, we regain the usual Poisson-Boltzmann equation, 
V 2 tp(r) cx exp [— /3qip(r)]. 

Since the boundary conditions are given in terms of 
the electric field strength, 



E{\r\=a) = -^L and £(| r | = R n ) = 



(30) 



to perform the numerical calculations it is convenient to 
rewrite the equations in terms of this variable. We now 
take advantage of the spherical symmetry of the system 
to eliminate the angular dependence of the equations, 
that is, we replace r by r = \r\ in Eqs. ( p6| ) to j2S|). The 
electric field also has only a spherically symmetric radial 
component, so that E(r) = E(r)y. Integrating the Pois- 
son equation ( p9| ) over a sphere of radius r and using the 
divergence theorem, we obtain a relation for the electric 
field strength E(r), 
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Inserting (|26|) into (31) we obtain an integro-differential equation for the electric field, 
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where we have chosen the gauge in which i/>(r = a) = 0, 
and the density profile, p c {r), which is necessary to eval- 
uate /i C orr(?~), i s a l so written in terms of the electric field, 
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The integro-differential equation was solved iteratively to 
obtain the electric field E(r). The charge density p c (r) 
is then calculated using Eq. (|33|). Finally, the internal 
free energy of a n-cluster can be expressed in terms of 
the charge density and the electric field to be 
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while the internal partition function of a n-cluster is 
C„=exp (-(3F™ n )- 



0.20 



VII. THERMODYNAMIC PROPERTIES 

The total Helmholtz free energy of the polyelectrolyte 
solution is a sum of the entropic and the electrostatic 
contributions, 

f(Pf,{Pn}) = Z miX + / PC + / PP + / CC • (35) 

Minimization of the total Helmholtz free energy under 
the constraints of fixed number of polyions and counte- 
rions leads to the law of mass action, 

po + np f = p n , (36) 

where the chemical potential of a specie s is p s = 
—df/dp s . This results in a set of Z coupled nonlinear 
algebraic equations for the densities p ni whose form is 
suitable to the use of an iterative method. Starting from 
a uniform distribution of clusters {/?„}, we were able to 
solve the coupled system numerically. A sample of the 
distributions obtained is presented in Fig. 1. Two fea- 
tures are worth remarking. The counterion condensation 
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FIG. 1. Cluster-density distribution {p n } for Z — 1000, 
volume fraction (j> — ^ira 3 N p /V = 0.01, and various values of 
temperature. From left to the right, the values of the reduced 
temperature are T* — 100, 50, and 20. 

is more effective as the temperature decreases, and the 
width of the distribution is not very sensitive to the vari- 
ations in temperature. The pressure can be obtained as a 
Legendre transform of the Helmholtz free-energy density, 

P = f(Pf,{Pn}) +^2psPs ■ (37) 
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FIG. 2. Dependence of the dimensionless (total) pressure, 
(3pa 3 /Z, on the volume fraction tf> = ^ira 3 N p /V, for several 
values of Z(1000, 2000, 5000) and T* = 100. 
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FIG. 3. Dependence of the 



various contributions to 
" x on the vol- 
for Z = 1000, and 
(ideal-gas) contribu- 



the dimensionless (total) pressure, 
ume fraction <\> — ^Tva 3 N p /V 
T* = 100. Solid line: mixing 
tion; long-dashed line: polyion-counterion contribution; 
dot-dashed line: polyion-polyion contribution; dashed line: 
counterion-counterion contribution. 

In Fig. 2 we present the total pressure inside the poly- 
electrolyte solution, which is a monotonically increasing 
function of the density of polyions. In agreement with the 
current experimental and Monte Carlo results 12 1, no ev- 
idence of any phase transition is encountered. To allow 
for a better appreciation of the relative importance of 
all the terms, in Fig. 3 we present separately the various 
contributions to the total pressure. All the electrostatic 
terms give a negative contribution to the total pressure, 
what can be interpreted as a form of an induced effective 
attraction between all the particles. 



ACKNOWLEDGMENTS 

This work has been supported by the Brazilian agency 
CNPq (Conselho Nacional de Desenvolvimento Cientffico 
e Tecnologico) . 



[1] R. A. Marcus, J. Chem. Phys. 23, 1057 (1955); S. Alexan- 
der, P. M. Chaikin, P. Grant, G. J. Morales, P. Pincus, 
and D. Hone, J. Chem. Phys. 80, 5776 (1984); R. D. 
Groot, J. Chem. Phys. 95, 9191 (1991); T. Gisler, S. 
F. Schulz, M. Borkovec, H. Sticher, P. Schurtenberger, 
B. DAguanno, and R. Klein, J. Chem. Phys. 101, 9924 
(1994); V. Reus, L. Belloni, T. Zemb, N. Lutterbach, 
and H. Versmold, J. Phys. II France 7, 603 (1997); M. 
J. Stevens, M. L. Falk, and M. O. Robbins, J. Chem. 
Phys. 104, 5209 (1996); L. B. Bhuiyan and C. W. Out- 
hwaite, Mol. Phys. 87, 625 (1996); W. B. Russel and D. 
W. Benzing, J. Colloid Interface Sci. 83, 163 (1981); M. 
N. Tamashiro, Y. Levin, and M. C. Barbosa, to appear 
in J. Phys. II France. 

[2] H. Lowen, J.-P. Hansen, and P. A. Madden, J. Chem. 
Phys. 98, 3275 (1993). 

[3] M. E. Fisher and Y. Levin, Phys. Rev. Lett. 71, 3826 

(1993) ; Y. Levin and M. E. Fisher, Physica A 225, 164 
(1996); Y. Levin, X.-J. Li, and M. E. Fisher, Phys. Rev. 
Lett. 73, 2716 (1994); B. Guillot and Y. Guissani, Mol. 
Phys. 87, 37 (1996). See also G. R. Stell, K. C. Wu, and 
B. Larsen, Phys. Rev. Lett. 37, 1369 (1976); G. R. Stell, 
Phys. Rev. A 45, 7628 (1992). 

[4] N. Bjerrum, Kgl. Dan. Vidensk. Selsk. Mat.-Fys. Medd. 
7, 1 (1926); H. Falkenhagen and W. Ebeling, in Ionic In- 
teractions, S. Petrucci (ed.) (Academic Press, New York, 
1971), Vol. 1: but note slips in Eqs. (49) and (51). 

[5] (a) P. W. Debye and E. Huckel, Phys. Z. 24, 185 (1923); 
(b) A good review is available in D. A. McQuarrie, Sta- 
tistical Mechanics (Harper & Row, New York, 1976), 
Chap. 15. 

[6] A. Z. Panagiotopoulos, Fluid Phase Equil. 76, 97 (1992); 
92, 313 (1994); J. P. Valleau, J. Chem. Phys. 95, 584 
(1991). 

[7] (a) Y. Levin, M. C. Barbosa, and M. N. Tamashiro, Eu- 

rophys. Lett. 41, 123 (1998); (b) Y. Levin, Europhys. 

Lett. 34, 405 (1996); (c) Y. Levin and M. C. Barbosa, J. 

Phys. II France 7, 37 (1997). 
[8] L. Onsager, Chem. Rev. 13, 73 (1933). 
[9] B. V. Derjaguin and L. Landau, Acta Physicochimica 

(USSR) 14, 633 (1941); E. J. W. Verwey and J. Th. G. 

Overbeek, Theory of the Stability of Lyophobic Colloids 

(Elsevier, Amsterdam, 1948); M. Medina-Noyola and D. 

A. McQuarrie, J. Chem. Phys. 73, 6279 (1980); X.-J. 

Li, Y. Levin, and M. E. Fisher, Europhys. Lett. 26, 683 

(1994) ; M. E. Fisher, Y. Levin, and X.-J. Li, J. Chem. 
Phys. 101, 2273 (1994). 

[10] We shall leave untouched the still controversial question 
of the existence of short-ranged attractive interactions 



6 



between the two polyions [but see I. Sogami and N. Ise, J. 
Chem. Phys. 81, 6320 (1984), and Yan Levin (to be pub- 
lished)] . If this attraction exists, it will not affect strongly 
the thermodynamic properties of the solution, which are 
dominated by the counterions and their interactions with 
the polyions (see Fig. 3). 
[11] (a) S. Nordholm, Chem. Phys. Lett. 105, 302 (1984); (b) 



R. Penfold, S. Nordholm, B. Jonsson, and C. E. Wood- 
ward, J. Chem. Phys. 92, 1915 (1990). 
[12] V. Vlachy, C. H. Marshall, and A. D. J. Haymet, J. Am. 
Chem. Soc. Ill, 4160 (1989); T. Palberg and M. Wiirth, 
Phys. Rev. Lett. 72, 786 (1994); J. C. Crocker and D. G. 
Crier, Phys. Rev. Lett. 73, 352 (1994). 



